home *** CD-ROM | disk | FTP | other *** search
/ OpenGL Superbible (2nd Edition) / OpenGL SuperBible e2.iso / tools / GLUT-3.7 / LIB / GLE / intersect.h < prev    next >
Encoding:
C/C++ Source or Header  |  1998-08-12  |  13.7 KB  |  392 lines
  1. /*
  2.  * FUNCTION:
  3.  * This file contains a number of utilities useful to 3D graphics in
  4.  * general, and to the generation of tubing and extrusions in particular
  5.  * 
  6.  * HISTORY:
  7.  * Written by Linas Vepstas, August 1991
  8.  * Updated to correctly handle degenerate cases, Linas,  February 1993 
  9.  */
  10.  
  11. #include <math.h>
  12. #include "port.h"
  13. #include "vvector.h"
  14.  
  15. #define BACKWARDS_INTERSECT (2)
  16.  
  17. /* ========================================================== */
  18. /*
  19.  * the Degenerate_Tolerance token represents the greatest amount by
  20.  * which different scales in a graphics environment can differ before
  21.  * they should be considered "degenerate".   That is, when one vector is
  22.  * a million times longer than another, changces are that the second will
  23.  * be less than a pixel long, and therefore was probably meant to be
  24.  * degenerate (by the CAD package, etc.)  But what should this tolerance
  25.  * be?  At least 1 in onethousand (since screen sizes are 1K pixels), but
  26.  * les than 1 in 4 million (since this is the limit of single-precision
  27.  * floating point accuracy).  Of course, if double precision were used,
  28.  * then the tolerance could be increased.
  29.  * 
  30.  * Potentially, this naive assumption could cause problems if the CAD
  31.  * package attempts to zoom in on small details, and turns out, certain
  32.  * points should not hvae been degenerate.  The problem presented here
  33.  * is that the tolerance could run out before single-precision ran
  34.  * out, and so the CAD packages would perceive this as a "bug".
  35.  * One alternative is to fiddle around & try to tighten the tolerance.
  36.  * However, the right alternative is to code the graphics pipeline in
  37.  * double-precision (and tighten the tolerance).
  38.  *
  39.  * By the way, note that Degernate Tolerance is a "dimensionless"
  40.  * quantitiy -- it has no units -- it does not measure feet, inches,
  41.  * millimeters or pixels.  It is used only in the computations of ratios
  42.  * and relative lengths.
  43.  */
  44.  
  45. /* 
  46.  * Right now, the tolerance is set to 2 parts in a million, which
  47.  * corresponds to a 19-bit distinction of mantissas. Note that
  48.  * single-precsion numbers have 24 bit mantissas.
  49.  */
  50.  
  51. #define DEGENERATE_TOLERANCE   (0.000002)
  52.  
  53. /* ========================================================== */
  54. /* 
  55.  * The macro and subroutine INTERSECT are designed to compute the
  56.  * intersection of a line (defined by the points v1 and v2) and a plane
  57.  * (defined as plane which is normal to the vector n, and contains the
  58.  * point p).  Both return the point sect, which is the point of
  59.  * interesection.
  60.  *
  61.  * This MACRO attemps to be fairly robust by checking for a divide by
  62.  * zero.
  63.  */
  64.  
  65. /* ========================================================== */
  66. /*
  67.  * HACK ALERT
  68.  * The intersection parameter t has the nice property that if t>1,
  69.  * then the intersection is "in front of" p1, and if t<0, then the
  70.  * intersection is "behind" p2. Unfortunately, as the intersecting plane
  71.  * and the line become parallel, t wraps through infinity -- i.e. t can
  72.  * become so large that t becomes "greater than infinity" and comes back 
  73.  * as a negative number (i.e. winding number hopped by one unit).  We 
  74.  * have no way of detecting this situation without adding gazzillions 
  75.  * of lines of code of topological algebra to detect the winding number;
  76.  * and this would be incredibly difficult, and ruin performance.
  77.  * 
  78.  * Thus, we've installed a cheap hack for use by the "cut style" drawing
  79.  * routines. If t proves to be a large negative number (more negative
  80.  * than -5), then we assume that t was positive and wound through
  81.  * infinity.  This makes most cuts look good, without introducing bogus
  82.  * cuts at infinity.
  83.  */
  84. /* ========================================================== */
  85.  
  86. #define INTERSECT(sect,p,n,v1,v2)            \
  87. {                            \
  88.    gleDouble deno, numer, t, omt;            \
  89.                             \
  90.    deno = (v1[0] - v2[0]) * n[0];            \
  91.    deno += (v1[1] - v2[1]) * n[1];            \
  92.    deno += (v1[2] - v2[2]) * n[2];            \
  93.                                \
  94.    if (deno == 0.0) {                    \
  95.       VEC_COPY (n, v1);                    \
  96.       /* printf ("Intersect: Warning: line is coplanar with plane \n"); */ \
  97.    } else {                        \
  98.                             \
  99.       numer = (p[0] - v2[0]) * n[0];            \
  100.       numer += (p[1] - v2[1]) * n[1];            \
  101.       numer += (p[2] - v2[2]) * n[2];            \
  102.                             \
  103.       t = numer / deno;                    \
  104.       omt = 1.0 - t;                    \
  105.                             \
  106.       sect[0] = t * v1[0] + omt * v2[0];        \
  107.       sect[1] = t * v1[1] + omt * v2[1];        \
  108.       sect[2] = t * v1[2] + omt * v2[2];        \
  109.    }                            \
  110. }
  111.  
  112. /* ========================================================== */
  113. /* 
  114.  * The macro and subroutine BISECTING_PLANE compute a normal vector that
  115.  * describes the bisecting plane between three points (v1, v2 and v3).  
  116.  * This bisecting plane has the following properties:
  117.  * 1) it contains the point v2
  118.  * 2) the angle it makes with v21 == v2 - v1 is equal to the angle it 
  119.  *    makes with v32 == v3 - v2 
  120.  * 3) it is perpendicular to the plane defined by v1, v2, v3.
  121.  *
  122.  * Having input v1, v2, and v3, it returns a unit vector n.
  123.  *
  124.  * In some cases, the user may specify degenerate points, and still
  125.  * expect "reasonable" or "obvious" behaviour.  The "expected"
  126.  * behaviour for these degenerate cases is:
  127.  *
  128.  * 1) if v1 == v2 == v3, then return n=0
  129.  * 2) if v1 == v2, then return v32 (normalized).
  130.  * 3) if v2 == v3, then return v21 (normalized).
  131.  * 4) if v1, v2 and v3 co-linear, then return v21 (normalized).
  132.  *
  133.  * Mathematically, these special cases "make sense" -- we just have to
  134.  * code around potential divide-by-zero's in the code below.
  135.  */
  136.  
  137. /* ========================================================== */
  138.  
  139. #define BISECTING_PLANE(valid,n,v1,v2,v3)            \
  140. {                                \
  141.    double v21[3], v32[3];                    \
  142.    double len21, len32;                        \
  143.    double dot;                            \
  144.                                 \
  145.    VEC_DIFF (v21, v2, v1);                    \
  146.    VEC_DIFF (v32, v3, v2);                    \
  147.                                 \
  148.    VEC_LENGTH (len21, v21);                    \
  149.    VEC_LENGTH (len32, v32);                    \
  150.                                 \
  151.    if (len21 <= DEGENERATE_TOLERANCE * len32) {            \
  152.                                 \
  153.       if (len32 == 0.0) {                    \
  154.          /* all three points lie ontop of one-another */    \
  155.          VEC_ZERO (n);                        \
  156.          valid = FALSE;                        \
  157.       } else {                            \
  158.          /* return a normalized copy of v32 as bisector */    \
  159.          len32 = 1.0 / len32;                    \
  160.          VEC_SCALE (n, len32, v32);                \
  161.          valid = TRUE;                        \
  162.       }                                \
  163.                                 \
  164.    } else {                            \
  165.                                 \
  166.       valid = TRUE;                        \
  167.                                 \
  168.       if (len32 <= DEGENERATE_TOLERANCE * len21) {        \
  169.          /* return a normalized copy of v21 as bisector */    \
  170.          len21 = 1.0 / len21;                    \
  171.          VEC_SCALE (n, len21, v21);                \
  172.                                 \
  173.       } else {                            \
  174.                                 \
  175.          /* normalize v21 to be of unit length */        \
  176.          len21 = 1.0 / len21;                    \
  177.          VEC_SCALE (v21, len21, v21);                \
  178.                                 \
  179.          /* normalize v32 to be of unit length */        \
  180.          len32 = 1.0 / len32;                    \
  181.          VEC_SCALE (v32, len32, v32);                \
  182.                                 \
  183.          VEC_DOT_PRODUCT (dot, v32, v21);            \
  184.                                 \
  185.          /* if dot == 1 or -1, then points are colinear */    \
  186.          if ((dot >= (1.0-DEGENERATE_TOLERANCE)) ||         \
  187.              (dot <= (-1.0+DEGENERATE_TOLERANCE))) {        \
  188.             VEC_COPY (n, v21);                    \
  189.          } else {                        \
  190.                                    \
  191.             /* go do the full computation */             \
  192.             n[0] = dot * (v32[0] + v21[0]) - v32[0] - v21[0];    \
  193.             n[1] = dot * (v32[1] + v21[1]) - v32[1] - v21[1];    \
  194.             n[2] = dot * (v32[2] + v21[2]) - v32[2] - v21[2];    \
  195.                                 \
  196.             /* if above if-test's passed,             \
  197.              * n should NEVER be of zero length */        \
  198.             VEC_NORMALIZE (n);                    \
  199.          }                             \
  200.       }                             \
  201.    }                                 \
  202. }
  203.  
  204. /* ========================================================== */
  205. /*
  206.  * The block of code below is ifdef'd out, and is here for reference
  207.  * purposes only.  It performs the "mathematically right thing" for
  208.  * computing a bisecting plane, but is, unfortunately, subject ot noise
  209.  * in the presence of near degenerate points.  Since computer graphics,
  210.  * due to sloppy coding, laziness, or correctness, is filled with
  211.  * degenerate points, we can't really use this version.  The code above
  212.  * is far more appropriate for graphics.
  213.  */
  214.  
  215. #ifdef MATHEMATICALLY_EXACT_GRAPHICALLY_A_KILLER
  216. #define BISECTING_PLANE(n,v1,v2,v3)                \
  217. {                                \
  218.    double v21[3], v32[3];                    \
  219.    double len21, len32;                        \
  220.    double dot;                            \
  221.                                 \
  222.    VEC_DIFF (v21, v2, v1);                    \
  223.    VEC_DIFF (v32, v3, v2);                    \
  224.                                 \
  225.    VEC_LENGTH (len21, v21);                    \
  226.    VEC_LENGTH (len32, v32);                    \
  227.                                 \
  228.    if (len21 == 0.0) {                        \
  229.                                 \
  230.       if (len32 == 0.0) {                    \
  231.          /* all three points lie ontop of one-another */    \
  232.          VEC_ZERO (n);                        \
  233.          valid = FALSE;                        \
  234.       } else {                            \
  235.          /* return a normalized copy of v32 as bisector */    \
  236.          len32 = 1.0 / len32;                    \
  237.          VEC_SCALE (n, len32, v32);                \
  238.       }                                \
  239.                                 \
  240.    } else {                            \
  241.                                 \
  242.       /* normalize v21 to be of unit length */            \
  243.       len21 = 1.0 / len21;                    \
  244.       VEC_SCALE (v21, len21, v21);                \
  245.                                 \
  246.       if (len32 == 0.0) {                    \
  247.          /* return a normalized copy of v21 as bisector */    \
  248.          VEC_COPY (n, v21);                    \
  249.       } else {                            \
  250.                                 \
  251.          /* normalize v32 to be of unit length */        \
  252.          len32 = 1.0 / len32;                    \
  253.          VEC_SCALE (v32, len32, v32);                \
  254.                                 \
  255.          VEC_DOT_PRODUCT (dot, v32, v21);            \
  256.                                 \
  257.          /* if dot == 1 or -1, then points are colinear */    \
  258.          if ((dot == 1.0) || (dot == -1.0)) {            \
  259.             VEC_COPY (n, v21);                    \
  260.          } else {                        \
  261.                                    \
  262.             /* go do the full computation */             \
  263.             n[0] = dot * (v32[0] + v21[0]) - v32[0] - v21[0];    \
  264.             n[1] = dot * (v32[1] + v21[1]) - v32[1] - v21[1];    \
  265.             n[2] = dot * (v32[2] + v21[2]) - v32[2] - v21[2];    \
  266.                                 \
  267.             /* if above if-test's passed,             \
  268.              * n should NEVER be of zero length */        \
  269.             VEC_NORMALIZE (n);                    \
  270.          }                             \
  271.       }                             \
  272.    }                                 \
  273. }
  274. #endif
  275.  
  276. /* ========================================================== */
  277. /*
  278.  * This macro computes the plane perpendicular to the the plane
  279.  * defined by three points, and whose normal vector is givven as the
  280.  * difference between the two vectors ...
  281.  * 
  282.  * (See way below for the "math" model if you want to understand this.
  283.  * The comments about relative errors above apply here.)
  284.  */
  285.  
  286. #define CUTTING_PLANE(valid,n,v1,v2,v3)                \
  287. {                                \
  288.    double v21[3], v32[3];                    \
  289.    double len21, len32;                        \
  290.    double lendiff;                        \
  291.                                 \
  292.    VEC_DIFF (v21, v2, v1);                    \
  293.    VEC_DIFF (v32, v3, v2);                    \
  294.                                 \
  295.    VEC_LENGTH (len21, v21);                    \
  296.    VEC_LENGTH (len32, v32);                    \
  297.                                 \
  298.    if (len21 <= DEGENERATE_TOLERANCE * len32) {            \
  299.                                 \
  300.       if (len32 == 0.0) {                    \
  301.          /* all three points lie ontop of one-another */    \
  302.          VEC_ZERO (n);                        \
  303.          valid = FALSE;                        \
  304.       } else {                            \
  305.          /* return a normalized copy of v32 as cut-vector */    \
  306.          len32 = 1.0 / len32;                    \
  307.          VEC_SCALE (n, len32, v32);                \
  308.          valid = TRUE;                        \
  309.       }                                \
  310.                                 \
  311.    } else {                            \
  312.                                 \
  313.       valid = TRUE;                        \
  314.                                 \
  315.       if (len32 <= DEGENERATE_TOLERANCE * len21) {        \
  316.          /* return a normalized copy of v21 as cut vector */    \
  317.          len21 = 1.0 / len21;                    \
  318.          VEC_SCALE (n, len21, v21);                \
  319.       } else {                            \
  320.                                 \
  321.          /* normalize v21 to be of unit length */        \
  322.          len21 = 1.0 / len21;                    \
  323.          VEC_SCALE (v21, len21, v21);                \
  324.                                 \
  325.          /* normalize v32 to be of unit length */        \
  326.          len32 = 1.0 / len32;                    \
  327.          VEC_SCALE (v32, len32, v32);                \
  328.                                 \
  329.          VEC_DIFF (n, v21, v32);                \
  330.          VEC_LENGTH (lendiff, n);                \
  331.                                 \
  332.          /* if the perp vector is very small, then the two     \
  333.           * vectors are darn near collinear, and the cut     \
  334.           * vector is probably poorly defined. */        \
  335.          if (lendiff < DEGENERATE_TOLERANCE) {            \
  336.             VEC_ZERO (n);                    \
  337.             valid = FALSE;                    \
  338.          } else {                        \
  339.             lendiff = 1.0 / lendiff;                \
  340.             VEC_SCALE (n, lendiff, n);                \
  341.          }                             \
  342.       }                             \
  343.    }                                 \
  344. }
  345.  
  346. /* ========================================================== */
  347.  
  348. #ifdef MATHEMATICALLY_EXACT_GRAPHICALLY_A_KILLER
  349. #define CUTTING_PLANE(n,v1,v2,v3)        \
  350. {                        \
  351.    double v21[3], v32[3];            \
  352.                         \
  353.    VEC_DIFF (v21, v2, v1);            \
  354.    VEC_DIFF (v32, v3, v2);            \
  355.                         \
  356.    VEC_NORMALIZE (v21);                \
  357.    VEC_NORMALIZE (v32);                \
  358.                         \
  359.    VEC_DIFF (n, v21, v32);            \
  360.    VEC_NORMALIZE (n);                \
  361. }
  362. #endif
  363.  
  364.  
  365. /* ============================================================ */
  366. /* This macro is used in several places to cycle through a series of
  367.  * points to find the next non-degenerate point in a series */
  368.  
  369. #define FIND_NON_DEGENERATE_POINT(inext,npoints,len,diff,point_array)   \
  370. {                                                                       \
  371.    gleDouble slen;                            \
  372.    gleDouble summa[3];                            \
  373.                                        \
  374.    do {                                                                 \
  375.       /* get distance to next point */                                  \
  376.       VEC_DIFF (diff, point_array[inext+1], point_array[inext]);        \
  377.       VEC_LENGTH (len, diff);                                           \
  378.       VEC_SUM (summa, point_array[inext+1], point_array[inext]);        \
  379.       VEC_LENGTH (slen, summa);                                         \
  380.       slen *= DEGENERATE_TOLERANCE;                    \
  381.       inext ++;                                                         \
  382.    } while ((len <= slen) && (inext < npoints-1));                      \
  383. }
  384.  
  385. /* ========================================================== */
  386.  
  387. extern int bisecting_plane (gleDouble n[3],    /* returned */
  388.                       gleDouble v1[3],  /* input */
  389.                       gleDouble v2[3],  /* input */
  390.                       gleDouble v3[3]);  /* input */
  391.  
  392.